Differential Abundance Analysis (DAA): Fecal microbiome

Based on the alpha and beta diversity results, I am not expecting to find a huge difference between groups in terms of bacterial composition. Please note that many of these code are from the Microbiome Intensive Course (MIC Course, Duke University, 2023) and thus the work of Dr. Pixu Shi.

This is after talking with Dr. Shi about what analyses I want to run.

1 Set-Up

First, update the ANCOMBC package if necessary.

#if (!requireNamespace("BiocManager", quietly=TRUE))
#     install.packages("BiocManager")
#BiocManager::install("ANCOMBC")

2 Load libraries and directories

# rm(list = ls(all = TRUE)) # Unload all packages if necessary
library(phyloseq)
library(fs)
library(rstatix)
library(ggpubr)
library(microViz)
library(ANCOMBC)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(DESeq2)
library(parallel)
library(dplyr)
library(tidyr)
library(readr)
library(tibble)
library(stringr)
library(cowplot)

rm(list=ls()) # clear environment 

git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME 
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 234 samples ]
## sample_data() Sample Data:       [ 234 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]

fig.dir = file.path(git.dir, "figures")

map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)

3 Prepare phyloseq object

I will now aggregate at the genus level, as it is the smallest taxonomic rank that analyses will be run on. In other words, I do not plan on running differential abundance analysis at ASV level.

Let’s prepare the phyloseq object for differential abundance analysis (DAA). Since we are interested in 1) differentially abundant taxa and 2) correlation analysis, we will limit the taxa to only ones that appear in > 3 samples for each sample type.

To avoid running analyses on extremely rare taxa that may be sequencing artifacts, the threshold for minimum number of counts is set at 10. This is, of course, arbitrary. 10 was chosen as 1% of number of read threshold of 1000 above for pruning samples. Based on how I wrote the code, taxa must be present have > 10 reads or more in > 3 samples.

Furthermore, categorical variables should be re-leveled as factors.

# Subset to true fecal samples only and remove tree
ps  %>% subset_samples(Sample_type == "Fecal" & Type == "True Sample") -> ps_fecal
ps_fecal_notree <- phyloseq(otu_table(ps_fecal), tax_table(ps_fecal), sample_data(ps_fecal))

# Categorical variables should be factors 
ps_fecal_notree@sam_data$Genotype <- factor(ps_fecal_notree@sam_data$Genotype, levels = c("WT", "KO"))
ps_fecal_notree@sam_data$Sex <- factor(ps_fecal_notree@sam_data$Sex, levels = c("Female", "Male"))
ps_fecal_notree@sam_data$Intranasal_Treatment <- factor(ps_fecal_notree@sam_data$Intranasal_Treatment, levels = c("PBS", "HDM"))
ps_fecal_notree@sam_data$Surgery <- factor(ps_fecal_notree@sam_data$Surgery, levels = c("None", "Sham", "VSG"))
ps_fecal_notree@sam_data$Mouse %>% unique() -> mouse_levels
ps_fecal_notree@sam_data$Mouse <- factor(ps_fecal_notree@sam_data$Mouse, levels = mouse_levels)
ps_fecal_notree@sam_data$Diet <- factor(ps_fecal_notree@sam_data$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
ps_fecal_notree@sam_data$Group <- factor(ps_fecal_notree@sam_data$Group, levels = c("NA", "Control", "1", "2"))
ps_fecal_notree@sam_data$Week <- factor(ps_fecal_notree@sam_data$Week, levels = c("0", "10", "13"))
ps_fecal_notree

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]

# Aggregate at genus level
ps_g_notree_fecal <- tax_glom(ps_fecal_notree, "Genus")
taxa_names(ps_g_notree_fecal) <- tax_table(ps_g_notree_fecal)[, 'Genus']

ps_g_notree_fecal

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 314 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 314 taxa by 6 taxonomic ranks ]

4 Getting ready: prune data

We will prune the phyloseq object at each phylogenetic rank at and above genus. Since we plan on running correlation analysis after differential abundance, let’s test for taxa that are present in at least 20% of all samples (arbitrary) with greater than 10 reads in each sample (also arbitrary. It’s also 1% of 1000 reads, which was the threshold used to remove samples from analysis).

x` Genus

# Pruning: set threshold number 
threshold_num = 0.20 * nrow(sample_data(ps_fecal))
nreads_prune = 10

# Prune taxa at each phylogenetic rank 
ps_g_notree_fecal %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_genus.prune
ps_genus.prune

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 56 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 56 taxa by 6 taxonomic ranks ]

ps_g_notree_fecal %>% transform_sample_counts(function(x) x/sum(x)) -> ps_genus.ts
ps_genus.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 314 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 314 taxa by 6 taxonomic ranks ]

subset_taxa(ps_genus.ts, ps_genus.ts@tax_table[, "Genus"] %in% ps_genus.prune@tax_table[, "Genus"]) -> ps_genus.prune.ts
ps_genus.prune.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 56 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 56 taxa by 6 taxonomic ranks ]

rowSums(ps_genus.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> genus_prune_RelAb

summary(genus_prune_RelAb$RelAb)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   67.11   99.32   99.65   99.25   99.85  100.00

4.1 Family

ps_fecal_notree %>% tax_glom("Family") -> ps_Family  
taxa_names(ps_Family) <- tax_table(ps_Family)[, 'Family']
ps_Family

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 169 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 169 taxa by 6 taxonomic ranks ]

ps_Family %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Family.ts
ps_Family.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 169 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 169 taxa by 6 taxonomic ranks ]

# Prune taxa at each phylogenetic rank 
ps_Family %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Family.prune
ps_Family.prune

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 28 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 28 taxa by 6 taxonomic ranks ]

subset_taxa(ps_Family.ts, ps_Family.ts@tax_table[, "Family"] %in% ps_Family.prune@tax_table[, "Family"]) -> ps_Family.prune.ts
ps_Family.prune.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 28 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 28 taxa by 6 taxonomic ranks ]

rowSums(ps_Family.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Family_prune_RelAb

summary(Family_prune_RelAb)

##      RelAb       
##  Min.   : 71.60  
##  1st Qu.: 99.68  
##  Median : 99.85  
##  Mean   : 99.54  
##  3rd Qu.: 99.94  
##  Max.   :100.00

4.2 Order

ps_fecal_notree %>% tax_glom("Order") -> ps_Order  
taxa_names(ps_Order) <- tax_table(ps_Order)[, 'Order']
ps_Order

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 109 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 109 taxa by 6 taxonomic ranks ]

ps_Order %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Order.ts
ps_Order.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 109 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 109 taxa by 6 taxonomic ranks ]

# Prune taxa at each phylogenetic rank 
ps_Order %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Order.prune
ps_Order.prune

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 18 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 18 taxa by 6 taxonomic ranks ]

subset_taxa(ps_Order.ts, ps_Order.ts@tax_table[, "Order"] %in% ps_Order.prune@tax_table[, "Order"]) -> ps_Order.prune.ts
ps_Order.prune.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 18 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 18 taxa by 6 taxonomic ranks ]

rowSums(ps_Order.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Order_prune_RelAb

summary(Order_prune_RelAb)

##      RelAb       
##  Min.   : 71.63  
##  1st Qu.: 99.76  
##  Median : 99.89  
##  Mean   : 99.56  
##  3rd Qu.: 99.96  
##  Max.   :100.00

4.3 Class

ps_fecal_notree %>% tax_glom("Class") -> ps_Class  
taxa_names(ps_Class) <- tax_table(ps_Class)[, 'Class']
ps_Class

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 42 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 42 taxa by 6 taxonomic ranks ]

ps_Class %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Class.ts
ps_Class.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 42 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 42 taxa by 6 taxonomic ranks ]

# Prune taxa at each phylogenetic rank 
ps_Class %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Class.prune
ps_Class.prune

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 12 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 12 taxa by 6 taxonomic ranks ]

subset_taxa(ps_Class.ts, ps_Class.ts@tax_table[, "Class"] %in% ps_Class.prune@tax_table[, "Class"]) -> ps_Class.prune.ts
ps_Class.prune.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 12 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 12 taxa by 6 taxonomic ranks ]

rowSums(ps_Class.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Class_prune_RelAb

summary(Class_prune_RelAb)

##      RelAb       
##  Min.   : 99.99  
##  1st Qu.:100.00  
##  Median :100.00  
##  Mean   :100.00  
##  3rd Qu.:100.00  
##  Max.   :100.00

4.4 Phylum

ps_fecal_notree %>% tax_glom("Phylum") -> ps_Phylum  
taxa_names(ps_Phylum) <- tax_table(ps_Phylum)[, 'Phylum']
ps_Phylum

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]

ps_Phylum %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Phylum.ts
ps_Phylum.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]

# Prune taxa at each phylogenetic rank 
ps_Phylum %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Phylum.prune
ps_Phylum.prune

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 6 taxonomic ranks ]

subset_taxa(ps_Phylum.ts, ps_Phylum.ts@tax_table[, "Phylum"] %in% ps_Phylum.prune@tax_table[, "Phylum"]) -> ps_Phylum.prune.ts
ps_Phylum.prune.ts

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 6 taxonomic ranks ]

rowSums(ps_Phylum.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Phylum_prune_RelAb

summary(Phylum_prune_RelAb)

##      RelAb    
##  Min.   :100  
##  1st Qu.:100  
##  Median :100  
##  Mean   :100  
##  3rd Qu.:100  
##  Max.   :100

5 Visualization

Let’s visualize at Phylum level to see how much information was kept:

plot_bar(ps_genus.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Genus", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Family.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Family", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Order.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Order", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Class.prune.ts, x = "Label", fill="Class") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Class, fill=Class), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Class", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Phylum.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Phylum", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

6 Subset metadata

meta.df %>% filter(Type == "True Sample" & Sample_type == "Fecal") %>% 
  filter(Label %in% phyloseq::sample_data(ps_g_notree_fecal)$Label) -> meta.fecal
meta.fecal

7 Q1: Fecal microbiome, HFD, wk 10 & 13, Surgery * Week

Let’s see how many mice we have:

meta.fecal %>% filter(Diet == "High_Fat_Diet" & Week == "13") -> meta.f.hfd
meta.f.hfd %>% 
  dplyr::count(Surgery, Genotype) %>% 
  arrange(n)

meta.f.hfd %>% 
  dplyr::count(Surgery) %>% 
  arrange(n)

nrow(meta.f.hfd)

## [1] 45

meta.f.hfd %>% distinct(Mouse) %>% nrow()

## [1] 45

7.1 Make TSE objects

ps_genus.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Genus
tse.fecal.Q1.Genus

## class: TreeSummarizedExperiment 
## dim: 56 45 
## metadata(0):
## assays(1): counts
## rownames(56): Clostridium sensu stricto 1 Anaerotruncus ... Romboutsia
##   Bifidobacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Family.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Family
tse.fecal.Q1.Family

## class: TreeSummarizedExperiment 
## dim: 28 45 
## metadata(0):
## assays(1): counts
## rownames(28): Peptococcaceae Clostridiaceae ... Anaerovoracaceae
##   Bifidobacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Order.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Order
tse.fecal.Q1.Order

## class: TreeSummarizedExperiment 
## dim: 18 45 
## metadata(0):
## assays(1): counts
## rownames(18): Peptococcales Clostridiales ...
##   Peptostreptococcales-Tissierellales Bifidobacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Class.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Class
tse.fecal.Q1.Class

## class: TreeSummarizedExperiment 
## dim: 12 45 
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
##   Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Phylum.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week == "13")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Phylum
tse.fecal.Q1.Phylum

## class: TreeSummarizedExperiment 
## dim: 10 45 
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
##   Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(45): SRR29054910 SRR29054915 ... SRR29055190 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

7.2 ANCOM2

7.2.1 Genus

set.seed(123) # Set seed for reproducibility 
n_cl_requested = 24

# Genus level
fecal.Q1.Genus <- ANCOMBC::ancom(data=tse.fecal.Q1.Genus, assay_name="counts", tax_level="Genus",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q1.Genus = fecal.Q1.Genus$res
res.fecal.Q1.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Genus.sig
tab.fecal.Q1.Genus.sig 

7.3 Family

set.seed(123) # Set seed for reproducibility 

# Family level
fecal.Q1.Family <- ANCOMBC::ancom(data=tse.fecal.Q1.Family, assay_name="counts", tax_level="Family",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q1.Family = fecal.Q1.Family$res
res.fecal.Q1.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Family.sig
tab.fecal.Q1.Family.sig 

7.3.1 Order

set.seed(123) # Set seed for reproducibility 

# Order level
fecal.Q1.Order <- ANCOMBC::ancom(data=tse.fecal.Q1.Order, assay_name="counts", tax_level="Order",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q1.Order = fecal.Q1.Order$res
res.fecal.Q1.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Order.sig
tab.fecal.Q1.Order.sig 

7.3.2 Class

set.seed(123) # Set seed for reproducibility 

# Class level
fecal.Q1.Class <- ANCOMBC::ancom(data=tse.fecal.Q1.Class, assay_name="counts", tax_level="Class",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q1.Class = fecal.Q1.Class$res
res.fecal.Q1.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Class.sig
tab.fecal.Q1.Class.sig 

7.3.3 Phylum

set.seed(123) # Set seed for reproducibility 

# Phylum level
fecal.Q1.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q1.Phylum, assay_name="counts", tax_level="Phylum",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q1.Phylum = fecal.Q1.Phylum$res
res.fecal.Q1.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Phylum.sig
tab.fecal.Q1.Phylum.sig 

7.3.4 Aggregate Results

tab.fecal.Q1.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q1.Family.sig$Tax_rank <- "Family"
tab.fecal.Q1.Order.sig$Tax_rank <- "Order"
tab.fecal.Q1.Class.sig$Tax_rank <- "Class"
tab.fecal.Q1.Phylum.sig$Tax_rank <- "Phylum"

rbind(tab.fecal.Q1.Genus.sig,
      tab.fecal.Q1.Family.sig,
      tab.fecal.Q1.Order.sig,
      tab.fecal.Q1.Class.sig,
      tab.fecal.Q1.Phylum.sig) -> tab.fecal.Q1

tab.fecal.Q1$Test <- "ANCOM II, fecal microbiome: all HFD, wk 10 & 13. Structural zeroes from Surgery, adjust with Genotype"

tab.fecal.Q1 

7.4 DESeq2

7.4.1 Genus

ds.fecal.Q1_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q1.Genus, ~ Surgery * Genotype)

## converting counts to integer mode

dds.fecal.Q1_Genus <- DESeq2::DESeq(ds.fecal.Q1_Genus)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

DESeq2::summary(dds.fecal.Q1_Genus)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1       580367       2515      33645         12845      12897.        5729.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    56          0           37.4

resultsNames(dds.fecal.Q1_Genus)

## [1] "Intercept"              "Surgery_Sham_vs_None"   "Surgery_VSG_vs_None"   
## [4] "Genotype_KO_vs_WT"      "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"

results(dds.fecal.Q1_Genus, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q1_Genus, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q1.Genus.Res
fecal.Q1.Genus.Res$Comparison <- c("Surgery: VSG over None")

results(dds.fecal.Q1_Genus, name  = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Genus, name  = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Sham surgery * KO genotype")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res

results(dds.fecal.Q1_Genus, name  = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("VSG * KO genotype")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res

fecal.Q1.Genus.Res$Compare_taxon <- "Genus"

fecal.Q1.Genus.Res

fecal.Q1.Genus.Res %>% filter(row %in% tab.fecal.Q1.Genus.sig)

7.4.2 Family

ds.fecal.Q1_Family <- DESeq2::DESeqDataSet(tse.fecal.Q1.Family, ~ Surgery * Genotype)

## converting counts to integer mode

dds.fecal.Q1_Family <- DESeq2::DESeq(ds.fecal.Q1_Family)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 2 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

DESeq2::summary(dds.fecal.Q1_Family)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1       721636       3576      43678         15928      16036.        7515.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    28          0           21.3

resultsNames(dds.fecal.Q1_Family)

## [1] "Intercept"              "Surgery_Sham_vs_None"   "Surgery_VSG_vs_None"   
## [4] "Genotype_KO_vs_WT"      "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"

results(dds.fecal.Q1_Family, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Family.Res
fecal.Q1.Family.Res$Comparison <- c("Surgery: Sham over None")

results(dds.fecal.Q1_Family, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Family, name  = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Family, name  = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Family, name  = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) 

fecal.Q1.Family.Res$Compare_taxon <- "Family"

fecal.Q1.Family.Res

fecal.Q1.Family.Res %>% filter(row %in% tab.fecal.Q1.Family.sig$taxon)

7.4.3 Order

ds.fecal.Q1_Order <- DESeq2::DESeqDataSet(tse.fecal.Q1.Order, ~ Surgery * Genotype)

## converting counts to integer mode

dds.fecal.Q1_Order <- DESeq2::DESeq(ds.fecal.Q1_Order)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

DESeq2::summary(dds.fecal.Q1_Order)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1       724646       3587      43755         15943      16103.        7577.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    18          0           13.9

resultsNames(dds.fecal.Q1_Order)

## [1] "Intercept"              "Surgery_Sham_vs_None"   "Surgery_VSG_vs_None"   
## [4] "Genotype_KO_vs_WT"      "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"

results(dds.fecal.Q1_Order, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Order.Res
fecal.Q1.Order.Res$Comparison <- c("Surgery:S Sham over none")

results(dds.fecal.Q1_Order, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")

results(dds.fecal.Q1_Order, name  = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Order, name  = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Order, name  = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05)

fecal.Q1.Order.Res$Compare_taxon <- "Order"

fecal.Q1.Order.Res

fecal.Q1.Order.Res %>% filter(row %in% tab.fecal.Q1.Order.sig$taxon)

7.4.4 Class

ds.fecal.Q1_Class <- DESeq2::DESeqDataSet(tse.fecal.Q1.Class, ~ Surgery * Genotype)

## converting counts to integer mode

dds.fecal.Q1_Class <- DESeq2::DESeq(ds.fecal.Q1_Class)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 2 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

DESeq2::summary(dds.fecal.Q1_Class)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1       729737       3587      43771         15951      16216.        7521.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    12          0           9.13

resultsNames(dds.fecal.Q1_Class)

## [1] "Intercept"              "Surgery_Sham_vs_None"   "Surgery_VSG_vs_None"   
## [4] "Genotype_KO_vs_WT"      "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"

results(dds.fecal.Q1_Class, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q1_Class, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q1.Class.Res
fecal.Q1.Class.Res$Comparison <- c("Surgery: VSG over None")

results(dds.fecal.Q1_Class, name  = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Class, name  = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Class, name  = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) 

fecal.Q1.Class.Res$Compare_taxon <- "Class"

fecal.Q1.Class.Res

tab.fecal.Q1.Class.sig

7.4.5 Phylum

ds.fecal.Q1_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q1.Phylum, ~ Surgery * Genotype)

## converting counts to integer mode

dds.fecal.Q1_Phylum <- DESeq2::DESeq(ds.fecal.Q1_Phylum)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

DESeq2::summary(dds.fecal.Q1_Phylum)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1       729772       3587      43771         15953      16217.        7521.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    10          0           7.84

resultsNames(dds.fecal.Q1_Phylum)

## [1] "Intercept"              "Surgery_Sham_vs_None"   "Surgery_VSG_vs_None"   
## [4] "Genotype_KO_vs_WT"      "SurgerySham.GenotypeKO" "SurgeryVSG.GenotypeKO"

results(dds.fecal.Q1_Phylum, tidy = TRUE, name = c("Surgery_Sham_vs_None")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Phylum.Res
fecal.Q1.Phylum.Res$Comparison <- c("Surgery: Sham over None")

results(dds.fecal.Q1_Phylum, name = c("Surgery_VSG_vs_None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res

results(dds.fecal.Q1_Phylum, name  = c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q1_Phylum, name  = c("SurgerySham.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Sham surgery * KO genotype")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res

results(dds.fecal.Q1_Phylum, name  = c("SurgeryVSG.GenotypeKO"), tidy = TRUE) %>%
  filter(padj < 0.05) 

fecal.Q1.Phylum.Res$Compare_taxon <- "Phylum"

fecal.Q1.Phylum.Res

tab.fecal.Q1.Phylum.sig

7.4.6 Aggregate results

rbind(fecal.Q1.Genus.Res, 
      fecal.Q1.Family.Res,
      fecal.Q1.Class.Res,
      fecal.Q1.Order.Res,
      fecal.Q1.Phylum.Res) -> fecal.Q1.deseq

fecal.Q1.deseq$Test <- "DESeq2, Fecal microbiome: all HFD, wk 10 + 13, Surgery * Genotype"

fecal.Q1.deseq

7.5 Intersection

tab.fecal.Q1 %>% filter(taxon %in% fecal.Q1.deseq$row) -> ancom.intersect.hfd
ancom.intersect.hfd

fecal.Q1.deseq %>% filter(row %in% tab.fecal.Q1$taxon) 

7.6 Make dataframe for graphing

psmelt(ps_genus.prune.ts) -> Genus.prune.melt 
psmelt(ps_Family.prune.ts) -> Family.prune.melt
psmelt(ps_Order.prune.ts) -> Order.prune.melt
psmelt(ps_Class.prune.ts) -> Class.prune.melt
psmelt(ps_Phylum.prune.ts) -> Phylum.prune.melt

bind_rows(
  Genus.prune.melt,
  Family.prune.melt,
  Order.prune.melt,
  Class.prune.melt,
  Phylum.prune.melt
)  -> fecal.RA.Full

fecal.RA.Full$Surgery <- factor(fecal.RA.Full$Surgery, levels = c("None", "Sham", "VSG"))

fecal.RA.Full %>% 
  dplyr::select(OTU:Label, Sample_type:Mouse, Genotype:Surgery, Week) -> fecal.RA.Full

head(fecal.RA.Full)

tail(fecal.RA.Full)

fecal.RA.Full %>% filter(Diet == "High_Fat_Diet" & Week != "0") -> fecal.graph.hfd
head(fecal.graph.hfd)

tail(fecal.graph.hfd)

8 Q2: no surgery group, all weeks, Diet * Genotype

Let’s see how many mice we have:

meta.fecal %>% filter(Surgery == "None") -> meta.f.nosurgery
meta.f.nosurgery %>% 
  dplyr::count(Week, Diet, Genotype) %>% 
  arrange(n)

meta.f.nosurgery %>% 
  dplyr::count(Week, Diet) %>% 
  arrange(n)

I will run Diet * Week + Genotype.

meta.f.nosurgery %>% distinct(Mouse) %>% nrow()

## [1] 58

8.1 Make TSE objects

ps_genus.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Genus
tse.fecal.Q2.Genus

## class: TreeSummarizedExperiment 
## dim: 56 120 
## metadata(0):
## assays(1): counts
## rownames(56): Clostridium sensu stricto 1 Anaerotruncus ... Romboutsia
##   Bifidobacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Family.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Family
tse.fecal.Q2.Family

## class: TreeSummarizedExperiment 
## dim: 28 120 
## metadata(0):
## assays(1): counts
## rownames(28): Peptococcaceae Clostridiaceae ... Anaerovoracaceae
##   Bifidobacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Order.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Order
tse.fecal.Q2.Order

## class: TreeSummarizedExperiment 
## dim: 18 120 
## metadata(0):
## assays(1): counts
## rownames(18): Peptococcales Clostridiales ...
##   Peptostreptococcales-Tissierellales Bifidobacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Class.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Class
tse.fecal.Q2.Class

## class: TreeSummarizedExperiment 
## dim: 12 120 
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
##   Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

ps_Phylum.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Phylum
tse.fecal.Q2.Phylum

## class: TreeSummarizedExperiment 
## dim: 10 120 
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
##   Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

8.2 ANCOM II

8.2.1 Genus

set.seed(123) # Set seed for reproducibility 

# Genus level
fecal.Q2.Genus <- ANCOMBC::ancom(data=tse.fecal.Q2.Genus, assay_name="counts", tax_level="Genus",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q2.Genus = fecal.Q2.Genus$res
res.fecal.Q2.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Genus.sig
tab.fecal.Q2.Genus.sig 

8.2.2 Family

set.seed(123) # Set seed for reproducibility 

# Family level
fecal.Q2.Family <- ANCOMBC::ancom(data=tse.fecal.Q2.Family, assay_name="counts", tax_level="Family",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q2.Family = fecal.Q2.Family$res
res.fecal.Q2.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Family.sig
tab.fecal.Q2.Family.sig 

8.2.3 Order

set.seed(123) # Set seed for reproducibility 

# Order level
fecal.Q2.Order <- ANCOMBC::ancom(data=tse.fecal.Q2.Order, assay_name="counts", tax_level="Order",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q2.Order = fecal.Q2.Order$res
res.fecal.Q2.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Order.sig
tab.fecal.Q2.Order.sig 

8.2.4 Class

set.seed(123) # Set seed for reproducibility 

# Class level
fecal.Q2.Class <- ANCOMBC::ancom(data=tse.fecal.Q2.Class, assay_name="counts", tax_level="Class",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q2.Class = fecal.Q2.Class$res
res.fecal.Q2.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Class.sig
tab.fecal.Q2.Class.sig 

8.2.5 Phylum

set.seed(123) # Set seed for reproducibility 

# Phylum level
fecal.Q2.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q2.Phylum, assay_name="counts", tax_level="Phylum",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 

## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!

res.fecal.Q2.Phylum = fecal.Q2.Phylum$res
res.fecal.Q2.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Phylum.sig
tab.fecal.Q2.Phylum.sig 

8.2.6 Aggregate Results

tab.fecal.Q2.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q2.Family.sig$Tax_rank <- "Family"
tab.fecal.Q2.Order.sig$Tax_rank <- "Order"
tab.fecal.Q2.Class.sig$Tax_rank <- "Class"
tab.fecal.Q2.Phylum.sig$Tax_rank <- "Phylum"

rbind(tab.fecal.Q2.Genus.sig,
      tab.fecal.Q2.Family.sig,
      tab.fecal.Q2.Order.sig,
      tab.fecal.Q2.Class.sig,
      tab.fecal.Q2.Phylum.sig) -> tab.fecal.Q2

tab.fecal.Q2$Test <- "ANCOM II, fecal microbiome: no surgery, main effect Diet (structural zeroes) and week + genotype"

tab.fecal.Q2 

8.3 DESeq2

8.3.1 Genus

ds.fecal.Q2_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q2.Genus, ~ Diet * Week + Genotype)

## converting counts to integer mode

dds.fecal.Q2_Genus <- DESeq2::DESeq(ds.fecal.Q2_Genus)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

DESeq2::summary(dds.fecal.Q2_Genus)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1231191       2007      45648         9976.      10260.        5536.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    56          0           40.0

resultsNames(dds.fecal.Q2_Genus)

## [1] "Intercept"                         "Diet_High_Fat_Diet_vs_Normal_Chow"
## [3] "Week_10_vs_0"                      "Week_13_vs_0"                     
## [5] "Genotype_KO_vs_WT"                 "DietHigh_Fat_Diet.Week10"         
## [7] "DietHigh_Fat_Diet.Week13"

results(dds.fecal.Q2_Genus, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q2_Genus, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Genus.Res
fecal.Q2.Genus.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

results(dds.fecal.Q2_Genus, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

results(dds.fecal.Q2_Genus, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

fecal.Q2.Genus.Res$Compare_taxon <- "Genus"

fecal.Q2.Genus.Res

8.3.2 Family

ds.fecal.Q2_Family <- DESeq2::DESeqDataSet(tse.fecal.Q2.Family, ~ Diet * Week + Genotype)

## converting counts to integer mode

dds.fecal.Q2_Family <- DESeq2::DESeq(ds.fecal.Q2_Family)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

DESeq2::summary(dds.fecal.Q2_Family)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1566638       2590      55908        12286.      13055.        6960.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    28          0           23.2

results(dds.fecal.Q2_Family, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q2_Family, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Family.Res
fecal.Q2.Family.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

fecal.Q2.Family.Res$Compare_taxon <- "Family"

fecal.Q2.Family.Res

8.3.3 Order

ds.fecal.Q2_Order <- DESeq2::DESeqDataSet(tse.fecal.Q2.Order, ~ Diet * Week + Genotype)

## converting counts to integer mode

dds.fecal.Q2_Order <- DESeq2::DESeq(ds.fecal.Q2_Order)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

DESeq2::summary(dds.fecal.Q2_Order)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1589033       2592      56999        12546.      13242.        7087.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    18          0           16.0

results(dds.fecal.Q2_Order, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q2_Order, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Order.Res
fecal.Q2.Order.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res

results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res

results(dds.fecal.Q2_Order, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q2_Order, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res

fecal.Q2.Order.Res$Compare_taxon <- "Order"

fecal.Q2.Order.Res

8.3.4 Class

ds.fecal.Q2_Class <- DESeq2::DESeqDataSet(tse.fecal.Q2.Class, ~ Diet * Week + Genotype)

## converting counts to integer mode

dds.fecal.Q2_Class <- DESeq2::DESeq(ds.fecal.Q2_Class)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

DESeq2::summary(dds.fecal.Q2_Class)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1594382       2595      57079         12562      13287.        7112.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    12          0           10.6

results(dds.fecal.Q2_Class, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q2_Class, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Class.Res
fecal.Q2.Class.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res

results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res

results(dds.fecal.Q2_Class, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q2_Class, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res

fecal.Q2.Class.Res$Compare_taxon <- "Class"

fecal.Q2.Class.Res

8.3.5 Phylum

ds.fecal.Q2_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q2.Phylum, ~ Diet * Week + Genotype)

## converting counts to integer mode

dds.fecal.Q2_Phylum <- DESeq2::DESeq(ds.fecal.Q2_Phylum)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

DESeq2::summary(dds.fecal.Q2_Phylum)

## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1594436       2595      57083         12562      13287.        7112.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    10          0           8.79

results(dds.fecal.Q2_Phylum, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 

results(dds.fecal.Q2_Phylum, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Phylum.Res
fecal.Q2.Phylum.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res

results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res

results(dds.fecal.Q2_Phylum, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) 

results(dds.fecal.Q2_Phylum, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res

fecal.Q2.Phylum.Res$Compare_taxon <- "Phylum"

fecal.Q2.Phylum.Res

8.3.6 Aggregate results

rbind(fecal.Q2.Genus.Res, 
      fecal.Q2.Family.Res,
      fecal.Q2.Class.Res,
      fecal.Q2.Order.Res,
      fecal.Q2.Phylum.Res) -> fecal.Q2.deseq

fecal.Q2.deseq$Test <- "DESeq2, Fecal microbiome: no surgery, ~ Diet * Week + Genotype"

fecal.Q2.deseq

8.4 Intersection

tab.fecal.Q2 %>% filter(taxon %in% fecal.Q2.deseq$row) -> ancom.intersect.nosurgery
ancom.intersect.nosurgery

ancom.intersect.nosurgery %>% arrange(Tax_rank, taxon)

fecal.Q2.deseq %>% filter(row %in% tab.fecal.Q2$taxon) %>%
  dplyr::select(row, Comparison, Compare_taxon, Test)  %>%
  group_by(row) %>%
  mutate(n.ind = paste0(1:n()))  %>%
  pivot_wider(
    names_from = "n.ind",
    values_from = Comparison
  )  -> fecal.Q2.deseq.intersect.wide

fecal.Q2.deseq.intersect.wide

8.5 Graph

fecal.RA.Full %>% filter(Surgery == "None") -> fecal.graph.nosurgery
head(fecal.graph.nosurgery)

tail(fecal.graph.nosurgery)

for (i in 1:length(ancom.intersect.nosurgery$taxon)) {
  fecal.graph.nosurgery %>% 
  filter(OTU == ancom.intersect.nosurgery$taxon[i]) %>%  
  ggplot(aes(x = as.character(Week), y = Abundance, color = Diet)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Relative Abundance", title = str_glue(ancom.intersect.nosurgery$taxon[i], " Relative Abundance, fecal"),
       x="Week") + 
  facet_wrap(~Diet) + 
  theme_bw(base_size = 14) + 
  theme(legend.position = "top") -> p
  print(p)
}

8.6 Graph: by diet

for (i in 1:length(ancom.intersect.nosurgery$taxon)) {
  fecal.graph.nosurgery %>% 
  filter(OTU == ancom.intersect.nosurgery$taxon[i]) %>%  
  ggplot(aes(x = Diet, y = Abundance, color = Diet)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Relative Abundance", title = str_glue(ancom.intersect.nosurgery$taxon[i], " Relative Abundance, fecal"),
       x="Diet") + 
  facet_wrap(~Week) + 
  theme_bw(base_size = 14) + 
  theme(legend.position = "top") + 
  stat_compare_means(method = "wilcox.test", aes(label = after_stat(p.signif)),
                     label.x = 1.5) -> p
  print(p)
}

## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `data.frame()`:
## ! arguments imply differing number of rows: 0, 1

## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `data.frame()`:
## ! arguments imply differing number of rows: 0, 1

9 Genotype?

DESeq2 provides more information on the results than ANCOM II does in the sense that you can specifically set the contrast. There were a number of taxa that were differentially abundant based on genotype:

fecal.Q2.deseq %>% filter(row %in% tab.fecal.Q2$taxon) %>% 
  filter(str_detect(Comparison, "Genotype") == TRUE) -> nosurgery.genotype

nosurgery.genotype

fecal.Q2.deseq.intersect.wide %>% filter(row %in% nosurgery.genotype$row)

It looks like these taxa were also affected by diet and week.

for (i in 1:length(nosurgery.genotype$row)) {
  fecal.graph.nosurgery %>% 
  filter(OTU == nosurgery.genotype$row[i]) %>%  
  ggplot(aes(x = Genotype, y = Abundance, color = Genotype)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Relative Abundance", title = str_glue(nosurgery.genotype$row[i], " Relative Abundance, fecal"),
       x="Diet") + 
  facet_grid(vars(Diet), vars(Week)) + 
  theme_bw(base_size = 14)-> p
  print(p)
}

fecal.graph.nosurgery %>% 
  filter(OTU %in% nosurgery.genotype$row) %>%
  group_by(OTU, Diet, Week) %>% 
  filter(sum(Abundance) > 0) %>% # remove groups that are all zeroes 
  wilcox_test(Abundance ~ Genotype) %>%
  adjust_pvalue(method = "holm") %>%
  add_significance()

NS after multiple testing.

10 Export csv files

fecal.Q1.deseq %>% filter(row %in% ancom.intersect.hfd$taxon) -> fecal.Q1.deseq.intersect
fecal.Q2.deseq %>% filter(row %in% ancom.intersect.nosurgery$taxon) -> fecal.Q2.deseq.intersect

write_csv(fecal.RA.Full, file.path(git.dir, "fecal_relative_abundance.csv"), append = FALSE, col_names = TRUE)
#write_csv(ancom.intersect.hfd, file.path(git.dir, "fecal_ancom_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(ancom.intersect.nosurgery, file.path(git.dir, "fecal_ancom_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)
#write_csv(fecal.Q1.deseq.intersect, file.path(git.dir, "fecal_deseq_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(fecal.Q2.deseq.intersect, file.path(git.dir, "fecal_deseq_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)

11 Reproducibility

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.1               stringr_1.5.0              
##  [3] tibble_3.2.1                readr_2.1.4                
##  [5] tidyr_1.3.0                 dplyr_1.1.3                
##  [7] DESeq2_1.40.2               SummarizedExperiment_1.30.2
##  [9] Biobase_2.60.0              MatrixGenerics_1.12.3      
## [11] matrixStats_1.1.0           GenomicRanges_1.52.1       
## [13] GenomeInfoDb_1.36.4         IRanges_2.34.1             
## [15] S4Vectors_0.38.2            BiocGenerics_0.46.0        
## [17] plotly_4.10.3               DT_0.29                    
## [19] ggrepel_0.9.4               microbiome_1.22.0          
## [21] colorBlindness_0.1.9        RColorBrewer_1.1-3         
## [23] pals_1.8                    ANCOMBC_2.2.2              
## [25] microViz_0.11.0             ggpubr_0.6.0               
## [27] ggplot2_3.4.4               rstatix_0.7.2              
## [29] fs_1.6.3                    phyloseq_1.44.0            
## 
## loaded via a namespace (and not attached):
##   [1] bitops_1.0-7                   DirichletMultinomial_1.42.0   
##   [3] httr_1.4.7                     doParallel_1.0.17             
##   [5] numDeriv_2016.8-1.1            tools_4.3.1                   
##   [7] doRNG_1.8.6                    backports_1.4.1               
##   [9] utf8_1.2.4                     R6_2.5.1                      
##  [11] vegan_2.6-4                    lazyeval_0.2.2                
##  [13] mgcv_1.9-0                     rhdf5filters_1.12.1           
##  [15] permute_0.9-7                  withr_2.5.1                   
##  [17] gridExtra_2.3                  cli_3.6.1                     
##  [19] sandwich_3.1-0                 labeling_0.4.3                
##  [21] sass_0.4.7                     mvtnorm_1.2-3                 
##  [23] proxy_0.4-27                   yulab.utils_0.1.0             
##  [25] foreign_0.8-85                 dichromat_2.0-0.1             
##  [27] scater_1.28.0                  decontam_1.20.0               
##  [29] maps_3.4.1                     readxl_1.4.3                  
##  [31] rstudioapi_0.15.0              RSQLite_2.3.3                 
##  [33] generics_0.1.3                 gridGraphics_0.5-1            
##  [35] vroom_1.6.4                    gtools_3.9.4                  
##  [37] car_3.1-2                      Matrix_1.6-1.1                
##  [39] biomformat_1.28.0              ggbeeswarm_0.7.2              
##  [41] fansi_1.0.5                    DescTools_0.99.52             
##  [43] DECIPHER_2.28.0                abind_1.4-5                   
##  [45] lifecycle_1.0.3                multcomp_1.4-25               
##  [47] yaml_2.3.7                     carData_3.0-5                 
##  [49] rhdf5_2.44.0                   Rtsne_0.16                    
##  [51] grid_4.3.1                     blob_1.2.4                    
##  [53] crayon_1.5.2                   lattice_0.22-5                
##  [55] beachmat_2.16.0                mapproj_1.2.11                
##  [57] pillar_1.9.0                   knitr_1.44                    
##  [59] boot_1.3-28.1                  gld_2.6.6                     
##  [61] codetools_0.2-19               glue_1.6.2                    
##  [63] data.table_1.14.8              MultiAssayExperiment_1.26.0   
##  [65] vctrs_0.6.4                    treeio_1.24.3                 
##  [67] Rdpack_2.6                     cellranger_1.1.0              
##  [69] gtable_0.3.4                   cachem_1.0.8                  
##  [71] xfun_0.40                      rbibutils_2.2.16              
##  [73] S4Arrays_1.2.1                 survival_3.5-7                
##  [75] SingleCellExperiment_1.22.0    iterators_1.0.14              
##  [77] gmp_0.7-2                      TH.data_1.1-2                 
##  [79] nlme_3.1-163                   bit64_4.0.5                   
##  [81] bslib_0.5.1                    irlba_2.3.5.1                 
##  [83] vipor_0.4.5                    rpart_4.1.21                  
##  [85] colorspace_2.1-0               DBI_1.1.3                     
##  [87] Hmisc_5.1-1                    nnet_7.3-19                   
##  [89] ade4_1.7-22                    Exact_3.2                     
##  [91] tidyselect_1.2.0               bit_4.0.5                     
##  [93] compiler_4.3.1                 htmlTable_2.4.2               
##  [95] BiocNeighbors_1.18.0           expm_0.999-7                  
##  [97] DelayedArray_0.26.7            checkmate_2.3.0               
##  [99] scales_1.2.1                   digest_0.6.33                 
## [101] minqa_1.2.6                    rmarkdown_2.25                
## [103] XVector_0.40.0                 htmltools_0.5.6.1             
## [105] pkgconfig_2.0.3                base64enc_0.1-3               
## [107] lme4_1.1-34                    sparseMatrixStats_1.12.2      
## [109] fastmap_1.1.1                  rlang_1.1.1                   
## [111] htmlwidgets_1.6.2              DelayedMatrixStats_1.22.6     
## [113] farver_2.1.1                   jquerylib_0.1.4               
## [115] zoo_1.8-12                     jsonlite_1.8.7                
## [117] energy_1.7-11                  BiocParallel_1.34.2           
## [119] BiocSingular_1.16.0            RCurl_1.98-1.13               
## [121] magrittr_2.0.3                 Formula_1.2-5                 
## [123] scuttle_1.10.3                 GenomeInfoDbData_1.2.10       
## [125] Rhdf5lib_1.22.1                munsell_0.5.0                 
## [127] Rcpp_1.0.11                    ape_5.7-1                     
## [129] viridis_0.6.4                  CVXR_1.0-11                   
## [131] stringi_1.7.12                 rootSolve_1.8.2.4             
## [133] zlibbioc_1.46.0                MASS_7.3-60                   
## [135] plyr_1.8.9                     lmom_3.0                      
## [137] Biostrings_2.68.1              splines_4.3.1                 
## [139] multtest_2.56.0                hms_1.1.3                     
## [141] locfit_1.5-9.8                 igraph_1.5.1                  
## [143] ggsignif_0.6.4                 rngtools_1.5.2                
## [145] reshape2_1.4.4                 ScaledMatrix_1.8.1            
## [147] evaluate_0.22                  tzdb_0.4.0                    
## [149] nloptr_2.0.3                   foreach_1.5.2                 
## [151] purrr_1.0.2                    rsvd_1.0.5                    
## [153] broom_1.0.5                    Rmpfr_0.9-3                   
## [155] e1071_1.7-13                   tidytree_0.4.5                
## [157] viridisLite_0.4.2              class_7.3-22                  
## [159] gsl_2.1-8                      lmerTest_3.1-3                
## [161] memoise_2.0.1                  beeswarm_0.4.0                
## [163] cluster_2.1.4                  TreeSummarizedExperiment_2.8.0
## [165] mia_1.8.0